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We consider the nonlinear transverse magnetic moment that arises in the Meissner state of 
superconductors with a strongly anisotropic order parameter. We compute this magnetic moment 
as a function of applied field and geometry, assuming d-wave pairing, for realistic samples, finite 
in all three dimensions, of high temperature superconducting materials. Return currents, shape 
effects, and the anisotropy of the penetration depth tensor are all included. We numerically solve 
' the nonlinear Maxwell-London equations for a finite system. Results are discussed in terms of the 

^\ I relevant parameters. The effect, which is a probe of the order parameter symmetry in the bulk, not 

0^ ■ just the surface, of the sample should be readily measurable if pairing is in a d-wave state. Failure 

to observe it would set a lower bound to the s-wave component. 

a ■ 
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I. INTRODUCTION 

^ ' The question of the pairii^ state of high temperature superconductors (HTSC's) continues to baffle researchers 
in the field. A recent reviewtl of many of the experimental results closes after over eighty pages with the tentative 
conclusion that the only state compatible with all experiments would have to exhibit two separate transitions, in 
\^ , contrast with the single transition invariably observed in HTSC's. Worse, evem different experiments performed on 
' the very same single crystal sample seem to lead to contradictory conclusions.H This situation makes it particularly 
\^ desirable to find good probes that may as unambiguously as possible discern the structure of the order parameter 
GN (OP) function. 

One such probe is afforded by the nonlinear Maxwell-London|-electrodynamics of exotic pairing states in supercon- 
ductors in the Meissner state. It was pointed out four years agoH in the context of a two-dimensional model (HTSC's 
are nearly all highly anisotropic, layered materials) that the presence of order parameter lines of nodes on a cylindrical 
FermLaurface would lead to observable nonlinear effects in the electrodynamic properties. These effects were shown 
(— ( to betrQ potentially capable of yielding a signature for the structure of the order parameter function, and therefore 
Q for the nodal structure of the gap itself. The most obvious of them is the presence of a transverse magnetic moment 
O m± perpendicular to the direction of an applied magnetic field H lying in the x — y {a — h) plane (the z axis is 
^ ' taken to be along the c crystallographic direction, perpendicular to the stack of planes in the lattice structure). This 
transverse magnetic moment has an angular dependence, as the sample is rotated about the z axis, reflecting directly 
the periodicity of the energy spectrum, that is, of the square of the order parameter, on the azimuthal angle cf). This 
signal would be particularly prominent if the order parameter has nodes or very deep minima as a function of angle. 
The transverse moment can be measured directly or inferred from the torque it produces. 

Electrodynamic effects are of particular interest because they probe the superconducting properties over a scale 
given by the field penetration depth, while many other experimental methods, such as tunneling and Josephson 
junctions, probe only over the scale of the correlation length, which in HTSC's is over two orders of magnitude 
smaller. Thus, electrodynamics probes the bulk properties of the supercoiiducting material, not merely the surface. 
Some of the apparent conflicts among experimental results may be resolveclQ if the order parameter symmetries in the 
bulk and at the surface differ. 

For these reasons, investigations of the transverse magnetic moment effect have continued. For the model of a "slab" 
shaped sample, infinite in the x — y directions and of fixed thickness, and a cylindrical Fermi surface, the equations 
involved in the nonlinear Maxwell-London electrodynamics have been solved numericallyu for the experimentally 
relevant ranges of applied magnetic field and temperature. At jT = 0, for the same geometry and assumptions, a 
perturbative, but accurate, analytic solution has been obtainedO. In Ref. || results were obtained for both a pure 
d-wave state and a njuxture of s and d waves. The amplitude of the transverse magnetization at low temperatures 
was found to increaseu as approximately the second power of the applied field, and, sinc&iOne is dealing with an effect 
arising from a finite penetration depth, as the surface area of the sample. The resultsua indicated that, for a pure 
d-wave state, mj_ should be detectable by a SQUID magnetometer in single crystal samples, even though the applied 
field is limited, by the need to keep the sample in the Meissner regime, to the field of first flux penetration. 
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The results of only one experiment seeking to find the transverse magnetizatiorB effect have been published!. No 
evidence of a transverse magnetic moment having the proper periodicity was obtained, but the nature of the data 
and of its analysis was such that the absence of evidence could not be conclusively taken to be evidence of absence. 
The resolution of the measurements was, over most of the field range studied, comparable to the signal measured. At 
the higher fields used the signal was below, but onlsi by a factor of two, the theoretically expected results for a pure 
d-wave state, based on the infinite slab assumptionsQ discussed in the above paragraph, scaled to the finite size of the 
sample. Although a subsequent analysisB of the field dependence of the experimental signaH was also supportive of 
the inference that no nodes were observed, a completely definitive statement about the existence of nodes could not 
be made because of theoretical and experimental uncertainties. 

Although this paper deals with the theoretical questions it is useful to briefly mention here the main difficulty 
associated with the experiments: except in the ideal world of theorists where all samples, even when finite, are spheres 
or other highly symmetric bodies, there always is in any superconductor, a tcansverse magnetic moment because of 
demagnetization effectsQ. For the rectangular samples used in the experimentEl these effects are much larger than the 
nonlinear contribution. Of course, their main periodicity, tt, is twice that associated (7r/2) with the d-wave nonlinear 
effects. Fourier analysis of the data, as done in Rcf. H] docs filter out the main spurious demagnetization siffnal, but 
its 7r/2 harmonic still will be confounded with the nonlinear signal and, although it can partly be sorted outn because 
of its different, linear, magnetic field dependence, the noise introduced in the original signal by the Fourier transforms 
and subtiactions takes a heavy toll. The work in this paper is motivated in large part by experiments now being 
plannedllj on samples having approximately a fiat disk shape that minimizes the demagnetization effect for an applied 
field lying in the a — b plane . It seems to us that reaching conclusions from experimental data taken on samples 
lacking rotational symmetry would be difficult. 

The theoretical uncertainties are related to the use in the work discussed above of two approximations: first, the 
use of an infinite slab geometry neglects the contribution to mj^ of return currents flowing along the sides of the 
sample, parallel to the z-axis. Furthermore, the penetration depth for those currents, Ac is much larger than that in 
the a-b plane, Xat- Since m_L increases with penetration depth, the question of how the larger Ac must be included 
in the analysis arises. The effect of Ac was included in the analysis of the datala but in a purely heuristic way, which 
assumed that its influence was rather strong. Although this assumption was relaxed in a subsequent reanalysisH, the 
whole question remains a major obstacle in reaching conclusions from experimental evidence. A flnitc Ac implies also 
abandoning the nation of a cylindrical Fermi surface. Further, the very neglect of finite size effects is suspect: they 
have been showrO to affect the analysis of penetration depth experiments in rather large samples. The uncertainties 
related to these complications in analyzing experimental data in the context of the geometric approximations of Ref. 
^ are about as large as the experimental uncertainties and it might be argued that, in combination, are about as large 
as the discrepancy between theory and experiment. 

In this paper we therefore undertake the examination of Maxwell-London electrodynamics, including nonlinear 
contributions, for finite samples. We include in the study the effects of Ac, which affect the electromagnetic response 
even at linear order in the field, and we specifically consider geometries relevant to the experimental systems currently 
been studied. Because of these realistic assumptions, not only is an analytic solution unattainable even at zero 
temperature, but the numerical work involved is quite considerable. Our results confirm that a transverse magnetic 
moment should be observable in obtainable single crystal samples, if the bulk pairing is in a d-wave state. The effect 
discussed in this paper is a fingerprint for the existence of an OP with d-wave symmetry. Its experimental observation 
in suitable samples, would constitute very hard to refute evidence for d-wave pairing in the bulk. Failure to observe 
this effect would at the very least put a lower bound on the existence of an s-wave component. 

In the next Section we present the equations solved and we discuss the methods we use. In Section III, we present 
our results and predictions for experimental outcomes, assuming a d-wave order parameter is present. Our conclusions 
are given in the last Section. 



II. METHODS 



A. Maxwell-London electrodynamics 

The equations of superconducting electrodynamics are in textbooksEO and need no rederivation. We will merely 
introduce here our notation, and briefly recall a couple of easy to overlook points. For the linear case, a particularly 
clear discussion of the complications that occur when the penetration depth tensor is anisotropic is in Chapter 3 of 
Ref.( §. 

Outside the sample the current is j = and therefore in the steady state the field H satisfies the Maxwell equations 
V • H = and V x H = 0. One can therefore write H in terms of a scalar potential <& satisfying the Laplace equation: 
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H = -V$ 



(2.1a) 



V2$ = (2.1b) 
In the interior of the sample one must consider three fields: H, j and the "superfluid velocity field" v defined as: 

where x is the phas e of the superconducting order parameter, A the vector potential, and e the proton charge. We 
set ?i = fcs = 1. Eq.(2_^) is equivalent to Eq. (5.99) in Ref. ||. The field v has dimensions of momentum, not velocity, 
but the factors associated with the effective mass and its anisotropics are more conveniently dealt with by placing 
them elsewhere. The three equations that one requires for these three fields are first, the second London equation 
obtained by taking the curl of ( ^.2[ ): 

V X V = -H (2.3) 

c 

second, Ampere's law: 

V X H = — j (2.4) 

c 

and finally, a constitutive equation relating j and v which is discussed in general in the next subsection. For the usual 
linear case this relation is of course simply: 

j = -epv (2.5) 



where p is the superfluid density tensor. The nonlinear contribution is discussed below (see Eq. ( 2.12 )). p is related 
to the penetration depth tensor A, whose components in the diagonal representation are the square of the London 
penetration depths, by the relation: 



When these tensors are proportional to the identity, then one can combine Eqns.( p.3[ ) and (2.4) and find that any one 
of the three fields considered satisfies the vector Hclmholtz equation. But this is not trueEl when p is anisotropic. One 
still has, however, the completely general equation: 

V X V X V - ^j(v) (2.7) 



which is valid whatever the relation j(v) might be. Eq. ( p.7[ ) will be the basic equation we will consider here. 

These equations must be solved with appropriate boundary conditions. These are: first, at infinity (that is, very 
far away from the finite sample) H must reduce to the applied field. Second, deep inside the sample all fields must 
vanish. Finally, at the interface H must be continuousEHlEj and the component of j normal to the interface must 
vanish. —. 

Once the currents j inside the sample are known, the magnetic moment can be obtained by integration:Lj 



m = 



|-ydrr xj(v(r)) (2., 



It is convenient to rewrite this equation in terms of surface integrals. Using Eq. (2.4) and formulas for vector 
calculual3 we derive: 

-^Isd'SU{n-r) + ^JJ^Snx^ (2.9) 

Integration is performed over the sample surface S, r is the position vector for a point on the surface S, and n is the 
unit normal pointing outwards. 

In the linear regime, and for the case which is experimentally relevant here where one has a small but finite value 
of the ratio X/d between a typical penetration depth A and a characteristic sample dimension d, the components 
of m can generically be written ia-the form m — mo(l — a{X/d) + 0{X/d)^). Examples for values of the positive 
constant a are given in textbookslljO. One has, then, a reduction in the magnetic moment due to current penetration 
in the material. The nonlinear effects may be conveniently viewed for our purposes as anisotropic, field-dependent 
corrections to the values of the a coefficients. 



3 



B. The supercurrent 



To obtain the supercurrent response one has to specify the full constitutive relation j(v) as discussed above. The 



linear contribution has been given in Eqns. (2j5) and (2.6). The anisotropy in A is related to that of the Fermi surface 
and can be related to the effective mass tensorcl. The anisotropy of the order parameter considered in this work would 
not, if the Fermi surface were isotropic, lead to any anisotropics in the penetration depthu, at linear order in the 
magnetic field. Thus we have: 

A=/^ (2.10) 

where m is the effective mass tensor and n the carrier density. j— . 

The nonlinear corrections to j(v) were first discussed by BardeenEJ in the context of an s-wave superconductor. In 
that case, at T = 0, the current-velocity relation is linear for velocities less than a critical velocity Vc defined below. 
Nonlinear corrections to the current-velocity relation arise from the thermal population of quasiparticles. These 
corrections are smalH (cubic in the ratia v /vc) and in the regime of larger flow velocities vortex nucleation may occur 
before these effects become importantEj. On the other hand, in unconventional superconductors, particularly those 
with nodes in the OP, nonlinear corrections are significantly larger and qualitatively different than in a superconductor 
with an s-wave excitation gap. The region on the Fermi surface (FS) close to the nodes in the gap provides enhanced 
nonlinear response due to the higher population of quasiparticles. There the energy required to break electron pairs 
will be reduced and will vanish along the direction of the lines of nodes in the gap. Thus nonlinear corrections in 
the presence of nodes in the OP are nonvanishing even at T = and the nonlinear behavior exhibits anisotropy with 
respect to the relative angle between applied magnetic field and the lines of nodes of the gap function. 

For the purpose of computing for a crystallographycally strongly anisotropic HTSC, in the geometry considered 
here, we are concerned with the angular dependence of the OP in the plane of the applied field, the a — h plane. In 
most of the numerical calculations performed in this work, we consider an order parameter of the pure d-wave form: 

A = Ao sin(2(/)) (2.11) 

where 4> is the azimuthal angle referred to a node and Aq the amplitude. The periodicity of m^, which equals that of 
the energy, is therefore tt/2. It would be possible, as we shall see, to generalize the calculations to other (f) dependence 
for the OP if necessary. Possible dependence of A on the c direction is very uncertain because of multilayering effectsEI, 
but we believe it is unlikely to be important. This follows from the above physical considerations, and from the brief 
discussion of results for alternative forms in the next Section and Appendix A. On the other hand, we must properly 
take into account the strong anisotropy of A in the c direction, which is related to that of the Fermi surface. 

The assumption of isotropy in the a — b place deserves further discussion, since anisotropics in the in-plane pene- 
tration depth in HTSC's are known to exist.Ej Anisotropy in Xab will le ad to a contribution to of periodicity tt 
(instead of 7r/2 for the nonlinear signal), and (see the paragraph below ( ^.9[ )) down by a factor of ~ X/d (the ratio 
of a typical penetration depth to a typical sample dimension) from the longitudinal, linear magnetic moment. This 
contribution will be linear, not quadratic, in the field. Hence, from an experimental point of view, it is in effect a small 
correction to the effective demagnetization factor which is, although to a small extent, present even in macroscopically 
symmetric samples. Thus, Fourier analysis of the experimental signal and examination of the field dependence of the 
harmonics would separate this from the nonlinear effect. As for the theoretical results discussed here it is sufficient 
to interpret our symbol Xab for the a — 6 plane penetration depth (Ac is that along the c-direction) as the geometric 
mean of Xa and Xb- Finally, if because of orthcirhombic anisotropy the nodes are not precisely at tt/2 angles, this 



can be accounted for by adding a small constantu to (2.11) and therefore the tt/2 Fourier component of the response 
would be little affected. |-.|-. 

For an order parameter with nodes, it is knowroO that the details of the FS shape are not important for the 
nonlinear properties we deal with here: only the nodes and their symmetry matter. Thus, our chief concern about the 
FS is to describe the anisotropy in A. For this purpose, we have used in our calculations an ellipsoidal (of revolution) 



Fermi surface characterized by effective masses niab (in the a — b plane) and rric, as in Eq. (Al). We introduce the ratio 
S = {rric/mab)- We define a speed in terms of the Fermi energy as = 2ef /niab which facilitates comparison 
with two-dimensional results. 

The general expressionQfl for j(v), including nonlinear terms, can be written as: 



J = -eNf j d's n(s)v/[(v/ ' v) + 2 f{E{S,) + • v)] (2.12) 

where N f is the total density of states at the Fermi Level, n{s) is the density of states at point s at the Fermi 



surface, normalized to unity, v/(s) is the s-dependent Fermi velocity, / the Fermi function, and E = \/ + |A(s)|^ 
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In general, the integrals in Eq. (2. IS) can be evaluated only numerically, as was done in 2-d in Ref. However it 
was shown there that at the temperatures (about two Kelvin) where the experiments are performed, one is very near 
the zero temperature limit. Therefore we confine ourselves in this work to this limit, where, within the assumptions 
discussed above, it is possible to derive an analytic expression for j(v), valid in the limit where S is large. Having 
an analytic expression for j(v) makes the subsequent substantial numerical work easier. In a magnetic field (hence 
V 7^ 0) even at T = there exists a region with |A(s)| + v/ • v < where it is possible to have a quasiparticle 
popul ation. One can then perform, as in Ref. |^ an approximate but accurate integration of the general formula Eq. 



( ^.12 ), to obtain the nonlinear corrections at zero temperature (details are given in Appendix A). The resulting j(v) 



relation is valid under the assumptions for the FS and the OP discussed above, and for any three dimensional strongly 
anisotropic superconductor (independent of sample geometry). Introducing cartesian axes x'-y' fixed to the crystal 
along the nodal directions the expression we obtain is: 

9tt 

jx',y' = -epabVx',v'{l " \'Vx' ,v'\) (2.13a) 

J, ^ -ep,.,(l - ^ ""''^f ) (2.13b) 
32'i;c \vx'\ + \ vy>\ 

where Vc = Aq/w/. The components of the superfluid density tensor are pab =^Nfv'j, pc = Pab/S, in the a — b plane 
and along the c crystallographic axis respectively. Even for YBCO reportedc3 measured values for i5 = rric/mab = 
i^c/^ab)'^ are quite large, between 15 and 50. Omitted terms in the above relation are those involving higher powers 
of v/vc <C 1 (for the regime of experimental interest in the Meissner state) or supressed by the small parameter S~^. 

C. Geometry 

In this subsection we discuss the sample geometry considered in the calculation of the supercurrent response in the 
Meissner state. The dependence of measured quantities such as the local magnetic field, actual current distribution 
and magnetic moment on sample shape is one of the issues of particular interest in this work, as explained in the 
Introduction. 

The experimentally relevant regime is that of small but finite ratio of penetration depths to typical sample dimen- 



sion. As mentioned below Eqn.(2.9), the geometric dependence of the effect under investigation in this regimeO is 
proportional to these ratios. Thus, the coefficient of this finite size term must be accurately calculated. A "locally 
flat" approximation (assuming a purely exponential decay for the fields away from the surface) is not sufficient for 
this purpose, as one can check even in the linear case from the exact solution for a sphereE ^. It is therefore necessary 



to solve the complete boundary value problem, that is, to find the solution to the equations (2J) and (2/7) outside and 
inside of the sample. This is in principle a numerically difficult undertaking, since a priori it involves solving a system 
of partial differential equations with nonlinear terms in the entire three dimensional space, but with currents being 
confined to a very small region where great precision is required. Although sophisticated variable grid methods could 
perhaps be found to take care of these complications, the relatively symmetric shape of the required experimental 
samples allows for a simpler approach. 

Single crystals of HTSC materials are typically flat, much thinner in the c-direction than in the a — b plane. The 
magnetic field, we recall, is applied parallel to this plane along the x direction. The geometry and coordinate system 
we consider are sketched in Fig. |l|. Because of the complications associated with demagnetization factors, as described 
in the Introduction, experinaejats must be performed on single crystals of a highly symmetric shape. This is achieved 
by laser cutting, or shaving23 the crystals so that their cross section in the direction perpendicular to the z axis is 
circular. Their shape is, therefore, roughly a disk, thinning towards the edges because of mechanical disintegration 
associated with the cutting. They have smooth edges and can therefore be described as ellipsoids of revolution. This 
geometry has also a considerable computational adva ntage: the form of the solution outside the sample in the limit 



A = is known exactly since the potential <& (see (2.1)) satisfies trivial Neumann boundary conditions at the interface 
and can be found by electrostatic analogy. The solution contains a single parameter which is simply related to mg, 
the value of the longitudinal magnetic moment in the zero penetration limit. When the penetration depth is finite, 
the longitudinal moment does change, but its correct value can be determined from the boundary conditions and the 
solution inside through an iteration process as described in the next Section. But, as important as these simplifications 
are, they only go so far: the fundamental equation (|]^) is not separable in spheroidal coordinates. 

We therefore consider a supercondutor in the shape of a flat ellipsoid of revolution (an oblate spheroid) with the 
axis of revolution along the z-axis. Its major and minor semiaxes are denoted by A and C respectively, and we have 
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A > C for actual samples. We take (see Fig. |T|) a coordinate system fixed to the direction of the magnetic field, 
with its z-axis parallel to the c crystalographic direction of the superconductor (and parallel to the C semiaxis of an 
ellipsoid). The field is applied along the x-axis, and we picture the experiment as being performed by rotating the 
crystal about the z-axis. As the crystal is rotated the axes x — y remain fixed in space, and should not be confused 
with the previously introduced x' — y' coordinates, affixed to the crystal structure. 

Even with the nonseparability of the equations, it is still convenient from the point of view of fuifilling the interface 
boundary conditions, to introduce oblate spheroidal coordinates a, /?, ip. In the definition we usecS, they are related 
to cartesian coordinates by the transformation: 

X = / cosh a sin /3 cos (y9 (2.14a) 
y = / cosh a sin /3 sin iy9 (2.14b) 

z = / sinh Q cos /? (2.14c) 

where < a < cxd, 0</3<7r,0<(y5< 27r, and / is a focal length scale factor. The surface of an ellipsoid corresponds 
to a given value of a. 

For an ellipsoid of revolution (about its z axis) in the presence of uniform magnetic field Ha applied aliwig the 
direction that we denote as the x axis the magnetic potential in the outside region for A = has the formcJ $ = 
—Hax{l + (7(sinha)), where the first term is the potential for the uniform applied field Hq and the gradient of the 
second term is the field generated by the superconducting ellipsoid. Writing out in detail the function g one has: 

$ = - Ha f Pl{isin\i a) Pl{cofi 13) cos if ^ AifQ\{isin\i a) Pl{cofi 13) cosip (2.15) 

where Pi and Q\ are associated Legendre functions of the first and second kind respectively. The parameter Ai is 
determined from the boundary conditions. It is proportional to m||, the longitudinal magnetic moment: 

mil = \fA^ (2.16) 



Its value at A = is: 



AM = , -ggsinhao ^2.17) 

1 + 1/cosh ao — sinhao arctan(l/ sinhao) 



where a = ao is the value at the surface of the ellipsoid. Eqns.(^T7) and ( 2.16| ) can be combined to yield the usual 



expression for mo, the zero penetration depth, purely longitudinal, magnetic moment, in terms of the ellipsoidal 
demagnetization factors .cJ At finite A the value for Ai is obtained from the iterative-procedure discussed in the next 
subsection. In the spherical limit sinha — > oo we also recover the standard result" ^rrl. 

III. RESULTS 
A. Numerical procedure 



Let us discuss now the iterative procedure we use to solve (2/7), while avoiding to have to numerically solve the 
field equations in all space. For clarity, let us focus first on the case where the nonlinear effects are neglected but the 
penetration depths are fi nite. This has two effects on the solution outside: first, the value of Ai (or of m\\) deviates 
from the value given by ( p7| ). Secondly, the field acquires higher order multipoles, i.e. the field potential acquires 
additional terms: 

$ = + AifQ\{ismha)Pl{cosl3) costp + ^ A„/Qi^(z sinh a)P,J (cos /3) cos(/3 (3.1) 

n>l 



where $a is the potential corresponding to the applied field (the first term on the right side of ( 2.15| )) and the 



gradient of the rest is the field created outside the sample by the current distribution in it. The terms with n > 1 
vanish at A = (and also in the spherical case if A is isotropic.) The fields generated by Ai are not purely dipolar, 
since they have ellipsoidal symmetry. However, the dipole moment is determined by this parameter only. At the 



6 



first step of the iteration we solve the linear ver sion of (2.7) for the current field, assuming that the outside field 
is given by its A = limit, i.e. the gradient of ( ^.15 ) with ( 2.17 ). Parameter counting shows that in order to do 



so we must give up one boundary condition, and accordingly we temporarily sacrifice th e co ntinuity of the "radial' 



field (d component). From the resulting current distribution we compute m through (2.9). This is of course not 
the same as the input-moment. We then replace this obtained value in the external potential (through (2.16)) and 
repeat the procedurfj until the moment generated by the computed currents equals the input value. The iteration 
is considerably simplified by observing that seiieral terms in ( |2.9| ) vanish explicitly when A = and by making use of 
the fact that the penetration depths are smallEj compared to the sample dimensions. Once the iteration is concluded, 
one finds that Ha is continuous except in a very narrow band near the equator corresponding to a symmet ry h igher 



than dipolar. This can then be eliminated by adding small values for higher order A„'s to the expansion (3.1), but 
symmetry considerations and examination of Eq. ( |2.9D show that these additions do not affect the already determined 
value of the sample magnetic moment. It is easy and very instructive to verify analytically that this procedure recovers 
the known resuluij for the A dependent magnetic moment of a sphere with an isotropic A. 

The same procedure is used with the nonlinear terms, with only two important diff erenc es: first, to the field outside 



one must add a transverse dipole, i.e. a contribution of the form of the last term in ( 2.15 ) rotated 90 



$i = Ai^fQl{ismha)P^{cosl3)smif (3.2) 

where ^f^Ai± = in±. Of course this term does not exist when the penetration depths vanish, since the nonlinear 
effects are absent unless the field can penetrate the sample. In principle one should also consider higher order multipole 



terms, as in (3.1) but we have found that any such terms are below the precision level of our numerical results. Both 



components of the mome nt m are determined through the vector relation (2.9). The other important difference is 
the obvious one of using ( ^.13 ) for the j(v) relation. There are also practical differences, however, since in the linear 



equations the variable (p can be separated out while in the full nonlinear case all three coordinates are coupled. 

In the actual solution of the equations we use a relaxation method. We proceed in two steps: first we solve the 
linear problem, which involves only two variables, since then the if dependence of all quantities can be determined 
analytically. The iterated solution for that problem is then used as the initial guess in the full three- variable nonlinear 
problem. 

We discretize the differential equations expressed in ellipsoidal (oblate spheroidal) coordinates, on a three dimens- 
sional grid. An obvious advantage of the ellipsoidal grid is that it simplifies consideration of boundary conditions at 
the surface of the ellipsoid, given by equation (|]^). Numerically, it is more intricate to consider boundary conditions 
that involve derivatives, and accuracy is increased if the grid points are also boundary points. 

The discretization procedure we have used has an estimated error quadratic in the spacing between the grid points. 
All quantities involved in equations ( ^.7^ have definite parity with respect to the exchange z — >— z and it is sufficient 
to solve them only for one half of the sample. Accuracy is predominantly governed by the spacing between grid 
points along the a direction. In the actual numerical solution (for half the sample), we consider ellipsoids of different 
shapes (different C/A) and sizes (different ^4). It is appropriate to increase Ua proportionally to C/A. Denoting 
by Ua the number of grid points along the a direction, the smallest used was 100, and the largest 800, spaced 
in the region of nonvanishing currents. The number of grid points along the /3 and ip directions was n/j = 50 and 
riip = 30 respectively. Increasing these numbers by a factor of two gave only effects below the numerical accuracy 
attained, which is about two significant figures, an error much smaller that the uncertainty arising from the imprecise 
knowledge of the experimental values of the input parameters. 

As one of the checks on the accuracy of the algorithm we use, it is instructive to consider the case of an isotropic 
spherical sample, (both the sample and the Fermi surface are spheres, 5 = 1) in the linear regime, where the analytical 
solution is known. To ensure that we tested the same algorithm, we treated the sphere as the the spherical limit 
of ellipsoidal coordinates, which corresponds to sinha — > oo. We used sinha =1000, equivalent to eccentricity 
e = a/1 — {C/A)^=0.001 with the grid given by Ua — 200, ri/3 — 50 and — 30 grid points. We solved equations 



(2.7) in a spherical shell of thickness lOA and studied both the magnetic moment and the current distributions. The 



accuracy for the current on any grid point corresponded t o fo ur significant figures for the regio n w here currents are 



important. The magnetic moment calculated both from (|2.8| ) and from the surface integrals (|2.9| ) agrees with the 
exact result, including the correct finite penetration depth correction. 

We have also checked that our results for the longitudinal magnetic moment extrapolated to the zero penetration 
depth limit agre e w ith the known analytic result for ellipsoids. We have also verified that the magnetic moment 
calculated from ( |2.8|) agrees from that found from the surface integrals (|2.9|) . The latter procedure is, however, much 
more convenient. 
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B. Numerical results and discussion 



In performing the calculations and describing the results, it is convenient to introduce dimensionless quantities. 
Because of the shape of the samples, we use Xab as the unit of length. We then define the dimensionless fields: V, J, 
and h: 

V = — , h=— , J = -— — J, (3.3) 

Vc Ho 47rAa6 

where we have introduced a characteristic magnetic field Hq as: 

wher e (/)o is the flux quantum and = Vf / ttAq is the in-plane superconducting coherence length. The definition of 
( |3.4[ ) involves precisely the same numerical factors as that used in Ref. ^. The requ ired equations are easily rewritten 



in terms of these quantities. The boundary conditions for the velocity field in (2/7) are now: 



V X V = h (3.5) 

where the right hand side is the external dimensionless field at the surface of the ellipsoid and from now on the 
derivatives are with respect to dimensionless length. The remain ing b oundary condition that there is no normal 
comp onent of current at the surface is readily obtained from Eqs. ( ^.13 ). The relation between j and v in equations 
( ^.13 ) then becomes: 



(V X V X V),.,y = -V,',y'{l - — \V^>,y'\) (3.6a) 

(iVxVxV)..-V.(l-|^|l|^) (3,6b) 

where we recall S = {Xc/Xab)'^ = rric/mab- Equations ( ^ ) are transformed to ellipsoidal coordinates, as defined above. 
Expressions for the superfluid density tensor and the differential operator V x V x V, in ellipsoidal coordinates are 
included in Appendix B. 

In Figures || and || we show some of the results for the currents. These figures illustrate some of the physics, as 
well as the quality of the numerics. In Fig ^ we show the current jz going along the z direction in the x — y plane 
as a function of distance from the surface of the sample starting at the point with cartesian coordinates (0,^,0). 
This current is overwhelmingly determined by the usual linear response to the field. One can clearly see that its 
decay as a function of depth from the surface is governed by the Ac penetration depth, and not by Xab- The next 
Figure illustrates the difference between the linear and nonlinear components of the current. It shows the current, 
again as a function of depth, starting at the center of the top of the sample, i.e. a point with cartesian coordinates 
(0,0, C). The component jy, normal to the applied field, is very predominantly "linear" and one can see that this 
time the decay is governed by the in-plane penetration depth, as expected from the geometry. The component jx 
along the applied field, on the other hand, arises exclusively from nonlinear effects: symmetry considerations show 
that it vanishes in the linear limit. Its overall scale is down by a large factor (basically the ratio of longitudinal and 
transverse moments). One sees that even though the overall decay of jx is determined by the scale Xab, its behavior 
is very far from exponential: it changes sign as a function of depth. This can be readily understood: at the positions 
plotted, close to the center of the the top of the ellipsoid, jx is approximately proportional (from Ampere's law) to 
the derivative with respect to z of the anomalous component of the field, Hy. This field component nearly vanishes at 
the surface (it would vanish for a slab) and decreases exponentially at depths larger than Xab- Hence, at this position 
in the sample, Hy has an extremum as a function of depth, and its derivative with respect to z must change sign at 
some point, as we find. This shows the delicate intricacy of the nonlinear current patterns inside the material. 

Before proceeding with the detailed discussion of the dependence of our results on the relevant physical quantitites 
(i.e. size, shape, applied field, and penetration depths) we illustrate their general scope by describing our prediction 
for the transverse magnetic moment of a possible HTSC superconducting sample, assumed to be in a pure d-wave 
pairing state. This is done in Fig. ^ We show there results for m± as a function of applied field. The quantity shown 
is the maximum value of m± as the crystal is rotated. It is assumed that the sample is an ellipsoid with A = 2 mm 
and C = 0.1 mm. Material parameters are taken to be = 20A, Xab = I8OOA, and Ac — 9000A. The characteristic 
field Ho would be about 5800 gauss. One can see that the magnetic moments involved are readily accessible to 
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measurement. Predictions for samples of other sizes, shapes and material parameters can be conveniently extracted 
from the information given below. 

We have obtained results for the transverse magnetic moment for a wide range of the experimentally accessible 
values of the appropriate dimcnsionless parameters, which as we shall see, can be taken to be Ha/Ho, Xc/Xab, '^ab/^, 
and C/A, the aspect ratio of the ellipsoid of revolution. We did not consider in our study the "thin film" situation 
where the sample is so small or so thin that its relevant dimensions are comparable to or smaller than the corresponding 
penetration depth. This case would be of no interest since the nonlinear effects then are vanishingly smallcl, and it is 
excluded from the analysis that follows. 

We begin our general discussion of the results by performing some dimensional analysis. The quantity Anm^ has 
dimensions of magnetic field times volume. The expression 

Q-'-^ (3.7) 

where V is the sample volume, is therefore dimcnsionless. Even though Q is suitable for some purposes, it is more 
convenient to analyze the dependence of the results on sample size in a different way. The reason is that 4:Trm± 
does not scale with the sample volume butQ as its surface area. The coefficient of proportionality between transverse 
moment and area depends on the sample shape and one can express this dependence through the aspect ratioEa A/C. 
Since the area of an ellipsoid of revolution is ttA^ times a function of A/C, it is easier for the purpose of giving results 
in a form more accessible to experimentalists, to scale explicitly results for ellipsoids of the same aspect ratio (i.e. the 
same sha,pe) by a factor of the "disk" area S = ttA^ . The third length that goes into the volume factor in the units of 
m± iaj'u a penetration depth. Since we are dealing with rather large A/C ratios, with currents predominantly in the 
a — b plane, it is advantageous for our purposes to take this length into accout by writing out a factor of Xab explicitly. 
Thus we put, as a first step: 

Anm^i^) = 7(^, ^,H)SXabfW (3-8) 

Here and hereafter we denote the magnitude of the applied field simply by H, rather than Ha, as there is no longer 
a possibility of confusion. It is not surprising in view of the_above dimensional analysis that we find 7 to be a rather 
weak function of its first two arguments and independentcj of S and of Xab /A. The angular dependence of m± is 
given by the function /(^) in terms of the angle -0 between the applied field and a node. We normalize f{ip) so that 
its value is unity at its maximum. 

Our results for the angular dependence are shown in Figure |^. The points shown are the values obtained from our 
numerical calculation. The error bars indicate the numerical uncertainty. The solid line represents the analytic result 
for the two dimensional calculations in the slab caseaij, normalized in the same way. We see from the Figure that 
the shape of f{ip) is, within numerical uncertainty, the same as for the flat case, where, with the same normalization, 
/(i/') = SV^sinip cosip{cosip — sin ■!/;), (for < V' < 7''/2). This is important, becauseO the 7r/2 Fourier coefficient 
of f{ip), as analytically calculated from the above expression, is unity to three significant figures. We can therefore 



identify here the coefficient of /(ip) in (3.8) with the Fourier amplitude m± as introduced earlier in the paper. The 
results shown in this Figure were obtained at H/Hq = 0.1, S = 16, and A/C — 19, but similar results are obtained in 
all other cases studied. 

As can already be seen in Figure 0, the coefficient 7 depends strongly on its third argument, the applied field. We 
have studied the field range < H/Hq < 0.2 at 0.05 intervals. We expect that the nonlinear is proportional to 
, as in the slab case at zero temperature. We therefore conclude that a very convenient way of writing our results 
for the amplitude is: 

4nm^^M{^,^)^HXabS fi^). (3.9) 

O Aab -no 

Equation ( ^j| ) implies explicitly that m±/H is a function of field only through the ratio H/Hq. We have verified 
that, as expected, M is independent of the field. This can be seen in Fig. |6[ where we plot the quantity G = 
10^4TTm± / {H XabS) vs H. We see that our results, represented by the symbols with error bars, are on a straight line, 
which is indicated by a best fit. Except for the field, the parameter values are the same as in the previous Figure. 
As in the case of the angular dependence, this qualitative result holds in all cases studied. The slope of the best fit 
straight line in plots such as that in this figure is used to extract M. 

The geometric aspect ratios we have considered in the nonlinear case range from A/C =7 up to 19 (that is from 
sinha =0.14433 to 0.05270). Comparison of the results with largest eccentricity to previous work on the slab case 
illustrates the effect of the return currents. This comparison cannot be made by taking C — > because we must have 
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C » Xab- The central p orti on of a flat ellipsoid can be identified as a "slab". It is simpler to make the comparison 



in terms of Q defined in (IT). It follows from our dimensional analysis that: 



Q = 7|^^/W (3.10) 

The corresponding result for the slab is of the same form, and the coefficient may be extracted from the T = results 
of Ref. m When all the jrelevant factors, such as the different definition of Hq^ the normalization of /(?A), and the 
setting of the parameterO (slope of the OP near a node) to = 2 are taken into account, one finds that the slab 
value is 7 = 0.056. We obtain 7 — 0.092 for our flattest ellipsoid dX 5 — 16, a number only weakly dependeat on 5. 
Hence, the presence of return currents enhances the nonlinear effects but, as hinted by the recent reanalysisci of the 
experimental data, the enhancement is about three times smaller, in typical situations, than the 5 dependent factor 
postulated in Ref. |[ 

We next show the dependence of our results on the remaining parameters, A/C and Ac/ Xab- As explained above, our 
results are most useful to experimentalists if given in terms of the scaled dimensionless amplitude AA. We therefore 
summarize the depe nden ce of this quantity on the mentioned length ratios in Table I. This table can be used, in 



conjuction with Eq. (3.9) to compute the theoretical predictions, when planning experiments or comparing data and 
theory. The rows are the relevant material parameter, given as the square of the penetration depth ratio, and the 
columns are the crystal shape. For crystals that are not quite ellipsoidal, one should choose a. C /A ratio in the table 
so that the surface to volume ratio of the ellipsoid agrees with that of the crystal. As the results vary slowly across 
the rows and columns of the Table, interpolation and reasonable extrapolation are obviously feasible. The value of 
Ai in the Table must be multiplied by the square of the applied field and by the cross sectional area of the sample in 
the a — 6 plane, then divided by Hq (which is again a material-dependent parameter), to obtain the expected value of 
the magnetic moment amplitude to be observed. The Table covers anisotropy parameter 5 values up to 50. 

The dependence of M. on its arguments is weak, less important than the unavoidable uncertainty in the values of 
the experimental input parameters such as those that go into e.g. the determination of Hq. The slow dependence of 
on ^ can be roughly understood as follows: the portion of the current loop which is along the c direction does of 
course contribute to the magnetic moment. However, provided, as in the Table, that A >> Ac > Xab the exact value 
of Ac does not matter very much, since (as seen in Fig. ^) the decay of the nonlinear current is governed primarily 
by Xab- The nonlinear moment is then limited by the smaller of the penetration depths. The slight downtrend with 5 
may be due to changes in the complicated current patterns discussed in connection with Fig. H. The slight decrease 
with decreasing eccentricity is likely to be due to the influence of the jz component (dominated by linear effects) 
in a thicker crystal. The Table indicates that flatter crystals are actually preferable to harder to obtain thick ones. 
Again, however, a simple explanation in intuitive terms is hard to come by, since these results must be related to the 
intricacies of the nonlinear current patterns. 

Finally we consider the sensitivity of our results to the unknown dependence of the order parameter on the crys- 
tallographic c di rection. In Appendix A, we have calculated the coefficients of the nonlinear terms in the current 



(as in Eq. ( ^.13 )) for the simplest form of the appropriateE3 3d OP. We have verified that an increase occurs in the 



quantity M. but it is negated by a compensating increase in iJo, where then one should replace ^0 by some average of 
in and out of plane correlation lengths. We conclude that the influence of any such dependence can be neglected as 
compared with the uncertainty in the material parameters. 



IV. CONCLUSIONS 



The detailed calculations for realistic samples performed here strongly confirm that the transverse magnetization 
effect should be readily observed in available crystals if the order parameter has nodes of the form expected for 
d-wave OP symmetry. From our results, see e.g. Fig. ^, it follows that the transverse magnetic moment generated 
in a typical sample of about 2 mm diameter and 0.1 mm thickness (the precise thickness not being very important) 
should be Airm^ « 10~^gauss cm'^. These numbers are above what can be readily measured by standard experimental 
techniques: the rather high uncertainty levels quoted in past experimentscl arise from the geometrical problems 
described in the Introduction and also from easy to overcome uncertainties in the angular positioning the sample. 
We have seen that the return currents in the z-direction produce an enhancement of the effect, although not nearly as 
large as that estimated in Ref. ^. The effect of the boundary conditions and the finite size on the result is somewhat 
intricate and, perhaps not surprisingly, we have not been able to find a clear physical picture in terms of some 
nonlinear generalization of the demagnetization factors. Although our results for the influence of sample geometry 
and anisotropy in penetration depth were not amenable to simple summing up either, the tabulated values of the 
quantity A4 should facilitate experimental design and interpretation. 
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As briefly indicated in Appendix A, (see below Eq.(All)) it is easy at T = to take into account different shapes of 



the d-wave order parameter function close to the nodes. Our calculations can also be extended to other symmetries of 
the OP, such as an s + d or s + zd state, or to take into account dependence of the OP in the c-direction. In the worst 
case, even at finite temperatures, the functional j(v) could be, numerically evaluated and used in the calculations, 
although possibly at considerable cost in computer time, by means of a lookup table and interpolation scheme. 

We have not considered in this work the effect of impurities. This is unnecessary as it was shown in Rcf. ^ that the 
impurity concentration required to noticeably decrease the nonlinear Meissner effect at low temperatures is such that 
it would clearly reduce the transition temperature of the material. This result will not depend on the geometry of the 
sample. On the other hand, we have considered here only the case where the order parameter has nodes, not just dips. 
The nonlinear effects are in fact extremely sensible to the presence of nodes and m± will decrease substantially, in finite 
as in infinite sampleaj if an, s- wave component eliminates the nodes. The magnitude of this decrease can be gauged 
from the infinite slab workQ since, again, it should not be excessively sensitive to return currents or sample shape. 
Therefore, a negative experimental result could be used to put a lower bound on the amount of s-wave component 
present. 

The methods discussed here can be extended to other sample shapes, to finite temperatures, and to the computation 
of other measurable effects that arise from the same nonlinear phenomena. 
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APPENDIX A: CURRENTS 



We derive here the expression for j(v) at T = starting from equation ( 2.12| ). Our procedurejis a generalization 



valid in the case of interest where the anisotropy S is large, of the two dimesional derivationQ. We consider an 
ellipsoidal Fermi surface: 

where we introduce Vf as in Section II. To establish our notation we briefly review the first, linear term in (^.12 ): 



j = -eNf J d^s n(s)v/(v/ • v) (A2) 

which is due to the condensate contribution to the current. For its evaluation one uses the convenient standard 
textbook method of rescaling kz by a factor of 6^^^^. The components of v/ are then given by: 

{vf)x' = Vf sm9 costj) (A3a) 

(vf)yi = Vf sin 6 sin (j) (A3b) 

{vf)z> = ^cose (A3c) 
Vo 

where (0, (f), {(/) measured from the x' axis) are the spherical angles of a vector with components kx,ky, kz/6^^'^. The 



scalar product in (A2) can then be written as: 



-Vf -v — Vf (v^i sin 9 cos </> + Vy' sin 9sin(/)-\ — — cos 9) (A4) 
The integration can be performed by replacing J cPs n{s) by J^^ dcj) d9 sin9 / {At:) which yields 
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jx'.y' — C-Pab^x' .y' (-^5) 
jz = -epcVz (A6) 

With pab = ^Nfvj:, pc = Pab/5, as quoted in Section II, (we recall that Nf is the total density of states, for both 
spins). Thus, as is well known, the effect of the ellipsoidal Fermi surfacS is to rescale pc by a factor of 1/5 relative 

to Pab- 

The nonlinear corrections are contained in the second term of ( 2.12| ) which is due to the quasiparticle backflow. We 
denote this term by jgp. At T=0, one can perform the f integration: 

jgp--2eiV/ / d^desine/{An) V/ e(-v/ • v - |A((^)|)a/(v/ • v)2 - \A{4>f (A7) 



where A(0) is given by ( 2.11 ). In the Meissner regime, \vfv\ <^ Aq and the contribution to the quasiparticles arises 
from narrow wedges along the nodal regions, approximately described by the azimuthal angle ^ i^/^c) ^1. A 
superfluid velocity v will give rise to backflow currents because of its components along nodal regions, separated by 
an azimuthal angle tt/2. In two dimensions, v can beO uniquely decomposed into two "jets" vix' and V2y' along two 
nodal directions. In the present three-dimensional case, where the anisotropy S is large, one can proceed in a similar 
way. We decompose v into the sum of two jets each directed along a nodal direction and tilted by the same angle lo 
with respect to the z-axis. Thus we write v = v^'x' + Vy>y' + = Vi + V2. with: 

vi — Vx'x' + cotuj\vx'\z (A8) 
V2 = Vy' y' + cot w I Wy' I z ( A9) 

cot UJ = r (AlO) 

\Vx'\ + \Vy'\ 

The decomposition thus specified is unique, and it ensures that the component is distributed along x'-z and y'-z 
planes proportionally to the corresponding projections of v along the x' , and y' directions. It is easy to see however, 
that any other decomposition of Vz would lead to the same results for the nonlinear currents derived below, except 
for higher order corrections in which shall be neglected in any case. 

The total phase space contributing to the quasiparticle part of the current can be obtained by considering the 



effects of Vi and V2 separately. Let us consider the effect of Vi in (A7). Quasiparticle excitations are allowed in the 
region described hy (j) < (j>c, where (v/ • vi)^ = (Aosin(/)c)^ and 

„ v'ivl, , „ 2 \vx'\ 1 T T . , . 

= ^L^(sin2 6'+ ^sin6icos6icottJ— + - cos^ 6i cot^ w (All 
4Ag V(5 Vx' 

where we have approximated the order parameter in the nodal region by A((/)) « Ao20. It would be easy to write 
instead A(0) « Ao/it/), with p, being a free parameter representing the slope of the OP function near the node, as was 
done in Ref. ||. This can be viewed as modifying the characteristic field Hq (Eq. ( ^.4[ ) by a factor of p/2. 

The integrals involved in the calculation of the nonlinear contribution to the relation j(v) due to the Vi are then 
of the form 



2Ao / deain^e / dcj) Vfyj4>l - 02 (A12) 



/2i = 2Ao/ desin9coa9— d0 w/V^c - 0^ (A13) 
Jo V S J -4,^ 

After integration over the angles {6, (j)) we get: 



StT^ 2 vl, ,^ . COt^ UJ , 



32 ' Vr io 



12 



Izl 



16 



(A15) 



The contribution to the quasiparicle current due to V2 , is calculated in a precisely similar manner and it gives analogous 
expressions for 7^/ and 1^2 ■ Comb ining the contributions due to Vi, V2 and omitting, as indicated above, subleading 
terms of order (1/(5) in Eq. (A14) we get: 



(A16) 



jqpz = eNfVf 



TT 1 

— (i'^'C0tcj + i;2, cotcj) 



(A17) 



In the above expression we substitute cotw given in (AlO) and recover the nonlinear part of Eq. (2. IS) 
If, as an alternative, we consider an OP function of the naive 3-d formtj cx (fc^ — ky): 

A{e,(j)) = Aisin^ (61) sin(20) 
then all integrals can still be done and one obtains instead: 



Jqp X' 
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(A18) 



(A19) 



TT 1 

jqp z = eNfv} — — {vl, cot uj + v^, cot w) 



'< 5 u, 



(A20) 



where is now defined in terms of Ai. The important nonlinear coefficient in ( A19) is a factor of 4/3 larger than 
that in AI6. The results for M increase by the same factor. However, in Eq. ( |3.4[ ) the quantity ^0 should be replaced 
by some average such as (CabCc)^^'^ which is smaller by about the same amount. 



APPENDIX B: QUANTITIES IN SPHEROIDAL COORDINATES 



The superfluid density tensor, given in the obvious cartesian coordinates by a diagonal tensor with components 
{Pabi Pab, Pc)j IS Converted to oblate spheroidal coordinates by performing the appropriate transformation. One obtains: 



P = Pab 



sinh^ Q sin^ /3+(5~^ cosh^ Q cos^ /3 (1— (5^ ^ ) sinh a cosh ct sin /3 cos /3 q 

sinh"^ a+cos^ P sinh"^ a-fcos^ /3 

(1 — ^~"^) sinh a cosh a sin /3 cos f3 sinh^ a sin^ /3+cosh'^ a cos^ /3 

sinh*^ a+cos^ /3 sinh"^ a+cos^ /3 

1 







To solve equations (2.7) the expression V x V x v should be transformed to oblate spheroidal coordinates: 
/^(V X V X v)q = aodfSi3Va + aidpVa + a2d^^Va + azVa + o-idapvp + a^davp + aedf3Vp 



(Bl) 



/^(V X V X v)^ bodaf3Va + hOaVa + b2dfjVa + h-^Va + h^daaVp + h^daVp + 6; 



(B2) 



/^(V X V X v)0 Poda4>Va + Pld^Va + P2dl34>Vp + Pid^Vfj + PidaaV^, + p^dppV^ 

+ PedaV^, + pvdpv^ + psv^, 



(B3) 



we recall that f is a focal length scale factor. The coefficients a^, bi,pi are given by (using the abreviations t — sinh a, 
u = cosha, s = sin/3, c = cos/3, w = {t'^ + c^y/^): 
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The right hand side of ( p/^ ) has to be transformed to spheroidal coordinates and added to the expressions above. 



^ J. Annett, N. Goldenfeld, and A. J. Leggett, to appear in Physical properties of High Temperature Superconductors, World 

Scientific, Singapore, (1996). 
^ J. Buan et al preprint. 

3 S.K. Yip and J.A. Sauls, Phys. Rev. Lett. 69, 2264, (1992). 
*B.P. Stojkovic and O.T. Vails, Phys. Rev. B51, 6049, (1995). 
^ D. Xu, S.K. Yip and J.A. Sauls, Phys Rev. B51, 16233, (1995). 
® See for example, S.R. Bahcall, unpublished. 

Magnetic moment is magnetization times volume. As both quantites are given in "emu" in gaussian units, much confusion 

between them occurs in the literature. Here we will give results in terms of 4irm, in gauss times volume, as it is in fact done 

in the experimental work cited here. 
* J. Buan, B.P. Stojkovic, N.E. IsraelofT, C.C. Huang, A.M. Goldman, and O.T. Vails, Phys. Rev. Lett. 72, 2632, (1994). 
® T.P. Orlando and K.A. Delin, Foundations of Applied Superconductivity, Addison- Wesley, Reading, (1991). 

A.M. Goldman, private communication. 
" Z.H. Lin, G.C. Spalding, A.M. Goldman, B.F. Bayman, and O.T. Vails, Europhysics Lett. 32, 573, (1995). 
^■^ P.G. DeGennes, Superconductivity of Metals and Alloys Addison- Wesley, reading, (1989). 

J.R. Reitz, F.J. Milford, and R.W. Christy, Foundations of Electromagnetic Theory, Addison- Wesley, Reading, (1989), page 

329. 

^"^ J. D. Jackson, Classical Electrodynamics, John Wiley, New York, (1975), page 181. 
J. D. Jackson, op. cit. inside cover. 

A. Fetter and J.D. Walecka, Theory of Many Particle Systems, McGraw-Hill, New York, (1971). 

F. London, Superfluids, v 1, John Wiley, New York, (1950). 
J. Bardeen, Rev. Mod. Phys. 34, 667, (1962). 

" D.N. Basov et al, Phys. Rev. Lett. 74, 598, (1995). K. Zhang et al, Phys. Rev. Lett. 73, 2484, (1994). 
See e.g. N.M. Plakida, High Temperature Superconductivity, Springer, Berlin, (1995). 

In the limit X ^ d the nonlinear effects disappear very quickly, as the sample becomes essentially normal. 
■^■^ A. Bhattacharya, D. Grupp, A.M. Goldman, and U. Welp, submitted to Applied Physics Letters. 

G. Arfken, Mathematical Methods for Physicists, Academic Press, New York, (1970), Chapter 2. This subject was dropped 
in later editions. 

N.N. Lebedev, LP. Skalskaya, and Y.S. Uflyand, Worked Problems in Applied Mathematics, Dover, New York, (1979). 

L.D. Landau and E.M. Lifshitz, Electrodynamics of Continuous Media, Pergamon, Oxford, (1960). 
■^^ In the "slab" case this procedure would yield immediately the correct result, since then there is no "radial" field and no 

boundary condition is sacrificed. 
■^'^ Additional mathematical details will be given elsewhere. 

One could use the eccentricity instead, but C/A is more immediately accessible. 



14 



The computational systems used, on which this scahng has been verified, range in size all the way up to the actual size of 
the physical samples. 

M. Sigrist and T.M. Rice, Zeitschrift fur Physik B (Condensed Matter) 68, 9, (1987) 

A. Bhattacharya private communication. 

C.H. Choi and P. Muzikar, Phys. Rev. B39, 11296, (1989). 



FIC. 1. Geometry considered in this paper. The sample is represented by an oblate ellipsoid of revolution. The left half of 
this figure shows the top view, and the right half a side view. The x, y and z directions are fixed in space. The field is applied 
along the x axis, as indicated, while is along the y axis. The sample is rotated about the z direction. The long and short 
semiaxes values are called A and C in the text, respectively. 



FIC. 2. The cartesian component of the current, jz, as a function of distance from the surface, beginning a the point with 
cartesian coordinates x = 0, y = A, z = 0. The quantity plotted is the current normalized to its value at the surface. D is 
distance in units of \ab- D = {A — y)/Xab- This current is dominated by linear effects in the field. One can see that its decay 
is governed by Ac, which is Ac = lAai, for the example plotted here, where we have used H/Hq — 0.1 and A/C = 7. 



FIC. 3. The linear and nonlinear components of the current as a function of distance from the surface beginning at the point 
X = y — 0, z = C. The dotted line is jy normalized to its value at the surface. It is predominantly linear and its decay is 
governed by Xab- The solid line represents the nonlinear component along the field direction, also normalized to its own, 
much smaller, surface value. Its nontrivial behavior is discussed in the text. D = (C — z)/Xab and all other parameters are as 
in the previous Figure. 



FIC. 4. Predictions of our theory for a typical HTSC single crystal sample. The size, shape, and material parameters are 
indicated in the text. The quantity plotted is the maximum value of 4-KTn±, (in units of 10~® gauss cm^) as the sample is 
rotated, as a function of applied field in gauss. 



FIC. 5. The angular dependence of the transverse magnetization, normalized to unity at its maximum. The points are 
our numerical results at several angles. The solid line represents the analytic result for the angular dependence in the purely 
two-dimensional case, normalized in the same way. Within error bars the angular dependence is the same, although the 
amplitude changes. 



FIC. 6. Dependence of the transverse mo men t amplitude on the applied field. The quantity plotted vs dimensionless field 
H/Ho, isG = 10'^ 4TTm±/{HXabS) (see Eq. (jsj).) The symbols are our numerical results and the dashed line the best linear 
fit. The linear dependence of G on H means that A4 is independent of field, and hence that m± is quadratic in H. 



TABLE I. The dimensionless quantity M, defined by Eq. ([j.9|), as a function of the material parameter 5 and of sample 
shape. 



S = (Xc/Xab)^ 


A/C=19 


A/C=10 


A/C=7 


16 


0.061 


0.060 


0.059 


25 


0.058 


0.057 


0.056 


36 


0.057 


0.054 


0.052 


50 


0.055 


0.051 


0.049 
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